#####################################################################
# Table 3, Figures 5-6
#####################################################################

#####################################################################
# Preliminaries
#####################################################################

rm(list = ls())

library(xtable)
library(scinference)
library(limSolve)
library(synthdid)
library(writexl)

#####################################################################
# Functions
#####################################################################

source("code/auxiliary scripts/common_functions.R")

#####################################################################
# Calibration DGPs
#####################################################################

source("code/auxiliary scripts/calibration_dgps.R")

#####################################################################
# Table 3
#####################################################################

# Setup

set.seed(12345)

nreps_sdid  <- 200
nreps_li    <- 5000
nreps_ttest <- 5000

# Simulations and Table 3

print("LaTeX output corresponding to Table 3, printed separately for each DGP")

combined_results_K3 <- combined_results_K4 <- matrix(NA,9,15)
  
for (DGP in 1:9){
  
  begin <- Sys.time()
  
  results_K3_cov <- results_K3_leng <- results_K3_bias <- matrix(NA,nreps_ttest,3)
  results_K4_cov <- results_K4_leng <- results_K4_bias <- matrix(NA,nreps_ttest,3)
  results_li_cov <- results_li_leng <- results_li_bias <- matrix(NA,nreps_li,1)
  results_sdid_cov <- results_sdid_leng <- results_sdid_bias <- matrix(NA,nreps_sdid,1)
  
  for (r in 1:nreps_ttest){

    res_K3_temp    <- sim_one_sample(DGP,T0=T0_app,T1=T1_app,J=J_app,K=3,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=FALSE,permtest = FALSE)
    res_K4_temp    <- sim_one_sample(DGP,T0=T0_app,T1=T1_app,J=J_app,K=4,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=FALSE,permtest = FALSE)
    
    results_K3_bias[r,]  <- res_K3_temp$bias_all
    results_K3_leng[r,]  <- res_K3_temp$leng_all
    results_K3_cov[r,]   <- res_K3_temp$cov_all
    
    results_K4_bias[r,]  <- res_K4_temp$bias_all
    results_K4_leng[r,]  <- res_K4_temp$leng_all
    results_K4_cov[r,]   <- res_K4_temp$cov_all
    
  }

  for (r in 1:nreps_li){
    
    res_li_temp <- sim_one_sample(DGP,T0=T0_app,T1=T1_app,J=J_app,K=3,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=TRUE,sdid=FALSE,permtest = FALSE)
    
    results_li_bias[r]  <- res_li_temp$bias_all
    results_li_leng[r]  <- res_li_temp$leng_all
    results_li_cov[r]   <- res_li_temp$cov_all
    
  }
  
  for (r in 1:nreps_sdid){
    
    res_sdid_temp <- sim_one_sample(DGP,T0=T0_app,T1=T1_app,J=J_app,K=3,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=TRUE,permtest = FALSE)

    results_sdid_bias[r]  <- res_sdid_temp$bias_all
    results_sdid_leng[r]  <- res_sdid_temp$leng_all
    results_sdid_cov[r]   <- res_sdid_temp$cov_all

  }
  
  cov_sdid  <- mean(results_sdid_cov)
  cov_li    <- mean(results_li_cov)
  
  bias_sdid <- mean(results_sdid_bias)
  bias_li   <- mean(results_li_bias)
  
  leng_sdid <- mean(results_sdid_leng)
  leng_li   <- mean(results_li_leng)
  
  cf_K3_bias <- colMeans(results_K3_bias)
  cf_K4_bias <- colMeans(results_K4_bias)
  
  cf_K3_cov <- colMeans(results_K3_cov)
  cf_K4_cov <- colMeans(results_K4_cov)
  
  cf_K3_leng <- colMeans(results_K3_leng)
  cf_K4_leng <- colMeans(results_K4_leng)
  
  overall_K3 <- c(10*cf_K3_bias[1:2],10*bias_li,10*bias_sdid,10*cf_K3_bias[3],cf_K3_cov[1:2],cov_li,cov_sdid,cf_K3_cov[3],cf_K3_leng[1:2],leng_li,leng_sdid,cf_K3_leng[3])
  overall_K4 <- c(10*cf_K4_bias[1:2],10*bias_li,10*bias_sdid,10*cf_K4_bias[3],cf_K4_cov[1:2],cov_li,cov_sdid,cf_K4_cov[3],cf_K4_leng[1:2],leng_li,leng_sdid,cf_K4_leng[3])
  
  combined_results_K3[DGP,] <- overall_K3
  combined_results_K4[DGP,] <- overall_K4
  
  print(DGP)
  print(Sys.time()-begin)
  print(xtable(cbind(c(3,4),rbind(overall_K3,overall_K4)),digits=c(0,0,rep(2,15))),include.rownames = F)

}

# Save table in .xlsx format

results_table3 <- rbind(sprintf("%.2f",combined_results_K3[1,]),sprintf("%.2f",combined_results_K4[1,]),
                        sprintf("%.2f",combined_results_K3[2,]),sprintf("%.2f",combined_results_K4[2,]),
                        sprintf("%.2f",combined_results_K3[3,]),sprintf("%.2f",combined_results_K4[3,]),
                        sprintf("%.2f",combined_results_K3[4,]),sprintf("%.2f",combined_results_K4[4,]),
                        sprintf("%.2f",combined_results_K3[5,]),sprintf("%.2f",combined_results_K4[5,]),
                        sprintf("%.2f",combined_results_K3[6,]),sprintf("%.2f",combined_results_K4[6,]),
                        sprintf("%.2f",combined_results_K3[7,]),sprintf("%.2f",combined_results_K4[7,]),
                        sprintf("%.2f",combined_results_K3[8,]),sprintf("%.2f",combined_results_K4[8,]),
                        sprintf("%.2f",combined_results_K3[9,]),sprintf("%.2f",combined_results_K4[9,]))

data_results_table3 <- as.data.frame(results_table3)
colnames(data_results_table3) <- NULL
write_xlsx(data_results_table3,"tables/table3.xlsx")

#####################################################################
# Figures 5-6
#####################################################################

# Setup

nreps_T1 <- 10000

T1s <- seq(10,16,2)

all_DGPs <- array(NA,dim = c(length(T1s),2,9))

# Simulations

for (DGP in 1:9){
  
  overall_cov_T1 <- overall_leng_T1 <- matrix(NA,length(T1s),1)
  
  for (t in 1:length(T1s)){
    
    res_temp_T1_cov <- res_temp_T1_leng <- matrix(NA,nreps_T1,1)
    
    for (r in 1:nreps_T1){
      
      res_temp_T1  <- sim_one_sample(DGP,T0=T0_app,T1=T1s[t],J=J_app,K=3,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=FALSE,permtest=FALSE)
      
      res_temp_T1_cov[r]   <- res_temp_T1$cov_all[1]
      res_temp_T1_leng[r]  <- res_temp_T1$leng_all[1]
      
    }
    
    overall_cov_T1[t]   <- mean(res_temp_T1_cov)
    overall_leng_T1[t]  <- mean(res_temp_T1_leng)
    
  }
  all_DGPs[,1,DGP] <- overall_cov_T1
  all_DGPs[,2,DGP] <- overall_leng_T1
  # print(DGP)
}

# Figures

for (DGP in 1:9){
  
  overall_cov_T1  <- all_DGPs[,1,DGP]
  overall_leng_T1 <- all_DGPs[,2,DGP]
  
  graphics.off()
  pdf(paste("graphics/small_T1_cov_",DGP,".pdf",sep=""),pointsize=16,width=8.0,height=6.0)
  plot(NULL, xlab=expression(T[1]), ylab="Coverage", main=paste("DGP",DGP), col="blue", pch=".", xlim = c(min(T1s),max(T1s)), ylim=c(0.5,1))
  lines(T1s,overall_cov_T1,col="black",lwd=5, lty=1,type="b",pch=1)
  dev.off()
  
  graphics.off()
  pdf(paste("graphics/small_T1_leng_",DGP,".pdf",sep=""),pointsize=16,width=8.0,height=6.0)
  plot(NULL, xlab=expression(T[1]), ylab="Average length", main=paste("DGP",DGP), col="blue", pch=".", xlim = c(min(T1s),max(T1s)), ylim=c(0,1))
  lines(T1s,overall_leng_T1,col="black",lwd=5, lty=1,type="b",pch=1)
  dev.off()
  
}
